<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of tt_rc</title>
  <meta name="keywords" content="tt_rc">
  <meta name="description" content="[Y]=TT_RC(D,N,ARR,ELEM_FUN,EPS,[OPTIONS])">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../index.html">Home</a> &gt;  <a href="#">tt2</a> &gt; <a href="index.html">cross</a> &gt; tt_rc.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../index.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for tt2/cross&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>tt_rc
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>[Y]=TT_RC(D,N,ARR,ELEM_FUN,EPS,[OPTIONS])</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>function [y]=tt_rc(d,n,elem_fun,eps,varargin) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">[Y]=TT_RC(D,N,ARR,ELEM_FUN,EPS,[OPTIONS])
The TT-Renormalization-Cross algorithm
input is: the dimension of the array d, the 
size vector n, the function that computes the prescribed element
elem_fum
Default parameters
nswp=10 (number of DMRG sweeps )
rmax=1000 (maximal rank of the solution)
x0=[];  (initial approximation)
verb = [0,1,2] (verbosity level)
vec=false (vectorization of the element function)</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="../../tt2/@qtt_tucker/diag.html" class="code" title="function [qt]=diag(qt)">diag</a>	Diagonal of a matrix or diagonal matrix from a vector in QTT-Tucker</li><li><a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>	Frobenius norm of the QTT-Tucker</li><li><a href="../../tt2/@tt_matrix/diag.html" class="code" title="function [tt]=diag(tm)">diag</a>	Extract the diagonal of the TT-matrix</li><li><a href="../../tt2/@tt_matrix/norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>	Matrix norm of the TT-matrix</li><li><a href="../../tt2/@tt_matrix/size.html" class="code" title="function [sz] = size(tt)">size</a>	Mode sizes of the TT-matrix</li><li><a href="../../tt2/@tt_tensor/diag.html" class="code" title="function [tm]=diag(tt)">diag</a>	Constructs diagonal TT-matrix from TT-tensor</li><li><a href="../../tt2/@tt_tensor/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>	Frobenius norm of the TT-tensor</li><li><a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>	Left and right orthogonalization of the TT-format</li><li><a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>	Reshape of the TT-tensor</li><li><a href="../../tt2/@tt_tensor/size.html" class="code" title="function [sz] = size(tt,dim)">size</a>	Mode sizes of the TT-tensor</li><li><a href="../../tt2/@tt_tensor/tt_tensor.html" class="code" title="function t = tt_tensor(varargin)">tt_tensor</a>	TT-tensor constructor</li><li><a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>	Maximal volume submatrix in an tall matrix</li><li><a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>	Truncation by absolution precision in Frobenius norm</li><li><a href="../../tt2/core/tt_ind2sub.html" class="code" title="function [ind] = tt_ind2sub(siz,ndx)">tt_ind2sub</a>	Correct conversion of an index to a multiindex</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="../../tt2/tests/test_cross.html" class="code" title="">test_cross</a>	Simple function to test the cross method</li><li><a href="../../tt2/tests/test_cross2.html" class="code" title="">test_cross2</a>	Simple function to test the cross method</li><li><a href="../../tt2/tests/test_cross3.html" class="code" title="">test_cross3</a>	</li></ul>
<!-- crossreference -->

<h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="#_sub1" class="code">function [super_core]=compute_supercore(i,elem_fun,d,n,ry,ind_left,ind_right,vec)</a></li><li><a href="#_sub2" class="code">function [iadd1,iadd2]=enlarge_cross(mat,i1,i2,eps)</a></li><li><a href="#_sub3" class="code">function [ind_left]=get_multi_left(i_left,r_left,ry)</a></li><li><a href="#_sub4" class="code">function [ind_right]=get_multi_right(i_right,r_right,ry)</a></li><li><a href="#_sub5" class="code">function [ind]=find_fiber_maximum(elem_fun,n,ind)</a></li></ul>

<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [y]=tt_rc(d,n,elem_fun,eps,varargin)</a>
0002 <span class="comment">%[Y]=TT_RC(D,N,ARR,ELEM_FUN,EPS,[OPTIONS])</span>
0003 <span class="comment">%The TT-Renormalization-Cross algorithm</span>
0004 <span class="comment">%input is: the dimension of the array d, the</span>
0005 <span class="comment">%size vector n, the function that computes the prescribed element</span>
0006 <span class="comment">%elem_fum</span>
0007 <span class="comment">%Default parameters</span>
0008 <span class="comment">%nswp=10 (number of DMRG sweeps )</span>
0009 <span class="comment">%rmax=1000 (maximal rank of the solution)</span>
0010 <span class="comment">%x0=[];  (initial approximation)</span>
0011 <span class="comment">%verb = [0,1,2] (verbosity level)</span>
0012 <span class="comment">%vec=false (vectorization of the element function)</span>
0013 
0014 <span class="keyword">if</span> ( numel(n) == 1 )
0015   n=n*ones(d,1);
0016 <span class="keyword">end</span>
0017 nswp=50;
0018 rmax=1000;
0019 x0=[];
0020 verb=2;
0021 vec=false; <span class="comment">%If there is a vectorization in options, i.e. if elem_fun</span>
0022 <span class="comment">%supports vectorized computations of many elements at once</span>
0023 change_dir_on=false;
0024 <span class="keyword">for</span> i=1:2:length(varargin)-1
0025     <span class="keyword">switch</span> lower(varargin{i})
0026         <span class="keyword">case</span> <span class="string">'nswp'</span>
0027             nswp=varargin{i+1};
0028         <span class="keyword">case</span> <span class="string">'rmax'</span>
0029             rmax=lower(varargin{i+1});
0030         <span class="keyword">case</span> <span class="string">'x0'</span>
0031             x0=varargin{i+1};
0032         <span class="keyword">case</span> <span class="string">'verb'</span>
0033             verb=varargin{i+1};
0034         <span class="keyword">case</span> <span class="string">'vec'</span>
0035             vec=varargin{i+1};
0036         <span class="keyword">case</span> <span class="string">'change_dir_on'</span>
0037             change_dir_on=varargin{i+1};
0038         <span class="keyword">otherwise</span>
0039             error(<span class="string">'Unrecognized option: %s\n'</span>,varargin{i});
0040     <span class="keyword">end</span>
0041 <span class="keyword">end</span>
0042 
0043 <span class="comment">%The initial setup is done.</span>
0044 
0045 <span class="comment">%We need initial index sets. Say, the vector of all ones (he-he)</span>
0046 <span class="comment">%i_left;</span>
0047 <span class="comment">%r_left;</span>
0048 <span class="comment">%i_right;</span>
0049 <span class="comment">%r_right;</span>
0050 i_left=cell(d,1);
0051 i_right=cell(d,1);
0052 r_left=cell(d,1);
0053 r_right=cell(d,1);
0054 
0055 <span class="keyword">if</span> ( isempty(x0) )
0056 ry=ones(d+1,1); <span class="comment">%Initial ranks</span>
0057 <span class="comment">%Find fiber maximum</span>
0058 ind=2*ones(d,1);
0059 <span class="comment">%ind=find_fiber_maximum(elem_fun,ind);</span>
0060 <span class="comment">%ind=randi([1,2],[d,1]);</span>
0061 ind=<a href="#_sub5" class="code" title="subfunction [ind]=find_fiber_maximum(elem_fun,n,ind)">find_fiber_maximum</a>(elem_fun,n,ind);
0062 <span class="keyword">for</span> i=1:d
0063   i_left{i}=ind(i);
0064   i_right{i}=ind(i);
0065   r_left{i}=1;
0066   r_right{i}=1;
0067 <span class="keyword">end</span>
0068 <span class="keyword">else</span>
0069   ry=x0.r;
0070   ps=x0.ps;
0071   cr=x0.core;
0072   rm=1;
0073   x0_old=x0;
0074   <span class="keyword">for</span> i=1:d-1
0075      cur_core=cr(ps(i):ps(i+1)-1);
0076      cur_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(i),n(i)*ry(i+1)]);
0077      cur_core=rm*cur_core;
0078      cur_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(i)*n(i),ry(i+1)]);
0079      [cur_core,rm]=<a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(cur_core,0);
0080      ind=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(cur_core); ind=ind';
0081      sbm=cur_core(ind,:);
0082      cur_core=cur_core / sbm; 
0083      <span class="comment">%q*inv(sbm)*sbm*rm</span>
0084      rm=sbm*rm;
0085      cr(ps(i):ps(i+1)-1)=cur_core(:);
0086      <span class="comment">%Convert ind to i_left r_left format</span>
0087      [radd,iadd]=ind2sub([ry(i);n(i)],ind);
0088      i_left{i}=iadd;
0089      r_left{i}=radd;
0090   <span class="keyword">end</span>
0091   cur_core=cr(ps(d):ps(d+1)-1);
0092   cur_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(d),n(d)*ry(d+1)]);
0093   cur_core=rm*cur_core;
0094   cr(ps(d):ps(d+1)-1)=cur_core(:);
0095   <span class="comment">%And backwards</span>
0096   rm=1;
0097   <span class="keyword">for</span> i=d:-1:2
0098      cur_core=cr(ps(i):ps(i+1)-1);
0099      cur_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(i)*n(i),ry(i+1)]);
0100      cur_core=cur_core*rm;
0101      cur_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(i),n(i)*ry(i+1)]);
0102      cur_core=cur_core.';
0103      [cur_core,rm]=<a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(cur_core,0);
0104      ind=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(cur_core);  ind=ind';
0105      sbm=cur_core(ind,:);
0106      cur_core=cur_core/sbm;
0107      rm=sbm*rm;
0108      cur_core=cur_core.';
0109      cr(ps(i):ps(i+1)-1)=cur_core;
0110      rm=rm.';
0111      [iadd,radd]=ind2sub([n(i);ry(i+1)],ind);
0112      i_right{i}=iadd;
0113      r_right{i}=radd;
0114   <span class="keyword">end</span>
0115   cur_core=cr(ps(1):ps(2)-1);
0116   cur_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(1)*n(1),ry(2)]);
0117   cur_core=cur_core*rm;
0118   cr(ps(1):ps(2)-1)=cur_core(:);
0119   x0.core=cr(:);
0120   <span class="comment">%keyboard;</span>
0121 <span class="keyword">end</span>
0122 
0123 <span class="comment">%Before we get the multiindex sets (helps to compute s-u-p-e-r-c-o-r-e-s)</span>
0124 
0125 ind_left=<a href="#_sub3" class="code" title="subfunction [ind_left]=get_multi_left(i_left,r_left,ry)">get_multi_left</a>(i_left,r_left,ry);
0126 ind_right=<a href="#_sub4" class="code" title="subfunction [ind_right]=get_multi_right(i_right,r_right,ry)">get_multi_right</a>(i_right,r_right,ry);
0127 
0128 <span class="comment">%Now we have to compute s-u-p-e-r-c-o-r-e-s (not cores!) using ind_left &amp;</span>
0129 <span class="comment">%ind_right</span>
0130 super_core=cell(d-1,1);
0131 <span class="keyword">for</span> i=1:d-1 <span class="comment">%There are d-1 s-u-p-e-r-c-o-r-e-s</span>
0132    <span class="comment">%Compute the index set for the i-th one</span>
0133    index_set=zeros(d,ry(i),n(i),n(i+1),ry(i+2));
0134    <span class="keyword">if</span> ( i == 1 )
0135        ileft=zeros(1,0);
0136    <span class="keyword">else</span>
0137       ileft=ind_left{i-1};    
0138    <span class="keyword">end</span>
0139    <span class="keyword">if</span> ( i == d - 1 )
0140       iright = zeros(1,0);
0141    <span class="keyword">else</span>
0142       iright=ind_right{i+2};
0143    <span class="keyword">end</span>
0144    <span class="keyword">for</span> s1=1:ry(i)
0145       <span class="keyword">for</span> s2=1:ry(i+2)
0146         <span class="keyword">for</span> i1=1:n(i)
0147            <span class="keyword">for</span> i2=1:n(i+1)
0148               index_set(:,s1,i1,i2,s2)=[ileft(s1,:),i1,i2,iright(s2,:)];
0149            <span class="keyword">end</span>
0150         <span class="keyword">end</span>
0151       <span class="keyword">end</span>
0152    <span class="keyword">end</span>
0153    M=ry(i)*n(i)*n(i+1)*ry(i+2);
0154    index_set=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(index_set,[d,M]);
0155    <span class="keyword">if</span> ( vec ) 
0156       super_core{i}=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(elem_fun(index_set),[ry(i),n(i),n(i+1),ry(i+2)]);
0157    <span class="keyword">else</span>
0158        cur_core=zeros(M,1);
0159       <span class="keyword">for</span> k=1:M
0160          cur_core(k)=elem_fun(index_set(:,k));
0161       <span class="keyword">end</span>
0162       super_core{i}=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(i),n(i),n(i+1),ry(i+2)]);
0163    <span class="keyword">end</span>
0164    
0165 <span class="keyword">end</span>
0166 <span class="comment">%Initialization is done.</span>
0167 <span class="comment">%Now to the main iteration</span>
0168 restart=true;
0169 <span class="keyword">while</span> ( restart )
0170     
0171 <span class="comment">%Through all the s-u-p-e-r-c-o-r-e-s</span>
0172 swp=1;
0173 dir=<span class="string">'lr'</span>;
0174 i=1;
0175 converged=false;
0176 rank_increase=false;
0177 <span class="keyword">while</span> ( swp &lt;= nswp &amp;&amp; ~converged)
0178    <span class="comment">%The method acts as follows.</span>
0179    <span class="comment">%1) Test if the current supercore is well approximated by the</span>
0180    <span class="comment">%current cross</span>
0181    <span class="comment">%2) If yes, do nothing.</span>
0182    <span class="comment">%3) If no, add the required indices to the cross and perform</span>
0183    <span class="comment">%modifications of the s-u-p-e-r-c-o-r-e-s</span>
0184    
0185    <span class="comment">%Nikrena ne soobrajau</span>
0186    <span class="comment">%Bilo: (i1,i2,i3,i4,i5)</span>
0187    <span class="comment">%Pofiksili (i1,i2), uvelichili r2. Rashirilos' mnojestvo (i2,i3,i4,i5)</span>
0188    <span class="comment">%(nikomu ne nujno)</span>
0189    <span class="comment">%V seredine: (i1,i2) (i3,i4,i5)+ - uvelichilos' mojectvo</span>
0190    <span class="comment">%Kogda uvelichili (i3,i4,i5) viros ry(3), uvelichilos 1 superyadro</span>
0191    cur_core=super_core{i}; 
0192    <span class="comment">%we have to convert the pair i_left{i} &amp; r_left{i} to the</span>
0193    <span class="comment">%index set for the matrix cur_core</span>
0194    
0195    cur_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(i)*n(i),n(i+1)*ry(i+2)]);
0196    <span class="comment">%i1=get_full_ind(i_left{i},r_left{i});</span>
0197    <span class="comment">%i2=get_full_ind(i_right{i+1},r_right{i+1});</span>
0198     <span class="comment">%ry(i)*n(i)</span>
0199     
0200     i1=r_left{i}+(i_left{i}-1)*ry(i); 
0201     i2=i_right{i+1}+(r_right{i+1}-1)*n(i+1);
0202  
0203   <span class="comment">%end</span>
0204   <span class="comment">%if ( swp == 3 &amp;&amp; i == 9 )</span>
0205   <span class="comment">%  keyboard;</span>
0206   <span class="comment">%end</span>
0207       <span class="comment">%If the cross is significantly &quot;bad&quot; then we can throw out certain</span>
0208       <span class="comment">%indices and redo, either recompute the full cross</span>
0209       [iadd1,iadd2]=<a href="#_sub2" class="code" title="subfunction [iadd1,iadd2]=enlarge_cross(mat,i1,i2,eps)">enlarge_cross</a>(cur_core,i1,i2,eps); 
0210       <span class="keyword">if</span> ( ~isempty(iadd1) ) <span class="comment">%There will be new elements in supercore{i+1} and supercore{i-1}!</span>
0211        rank_increase=true;
0212        <span class="comment">%Increase i_left{i} &amp; i_right{i+1}</span>
0213        [radd_left,iadd_left]=ind2sub([ry(i);n(i)],iadd1);
0214        [iadd_right,radd_right]=ind2sub([n(i+1);ry(i+2)],iadd2);
0215        i_left{i}=[i_left{i};iadd_left];
0216        r_left{i}=[r_left{i};radd_left];
0217        <span class="comment">%if ( i == 8 )</span>
0218        <span class="comment">%    keyboard</span>
0219        <span class="comment">%end</span>
0220        i_right{i+1}=[i_right{i+1};iadd_right];
0221        r_right{i+1}=[r_right{i+1};radd_right];
0222        rold=ry(i+1);
0223        ry(i+1)=ry(i+1)+numel(iadd1);
0224        <span class="comment">%Zaglushka</span>
0225        ind_left=<a href="#_sub3" class="code" title="subfunction [ind_left]=get_multi_left(i_left,r_left,ry)">get_multi_left</a>(i_left,r_left,ry);
0226        
0227        ind_right=<a href="#_sub4" class="code" title="subfunction [ind_right]=get_multi_right(i_right,r_right,ry)">get_multi_right</a>(i_right,r_right,ry);
0228        
0229        <span class="keyword">if</span> ( i &gt; 1 ) 
0230          <span class="comment">%Recompute supercore{i-1};</span>
0231          <span class="comment">%supercore{i-1}=zeros(ry(i-1),n(i-1),n(i),ry(i+1))</span>
0232          radd=numel(iadd1);
0233          rtmp=ry;  rtmp(i+1)=radd;
0234          <span class="comment">%tmp1=ind_left;</span>
0235          <span class="comment">%tmp1{i}(:,end-radd+1:end)=[];</span>
0236          tmp2=ind_right;
0237          tmp2{i+1}(1:rold,:)=[];
0238          core_add=<a href="#_sub1" class="code" title="subfunction [super_core]=compute_supercore(i,elem_fun,d,n,ry,ind_left,ind_right,vec)">compute_supercore</a>(i-1,elem_fun,d,n,rtmp,ind_left,tmp2,vec);
0239          new_core=zeros(ry(i-1),n(i-1),n(i),ry(i+1));
0240            new_core(:,:,:,1:rold)=super_core{i-1};
0241          new_core(:,:,:,rold+1:end)=core_add;
0242          super_core{i-1}=new_core;
0243          <span class="comment">%p1=compute_supercore(i-1,elem_fun,d,n,ry,ind_left,ind_right,vec);</span>
0244          <span class="comment">%keyboard;</span>
0245        <span class="keyword">end</span>
0246        <span class="keyword">if</span> ( i &lt; d - 1)
0247          radd=numel(iadd1);
0248          rtmp=ry;  rtmp(i+1)=radd;
0249          tmp1=ind_left;
0250          tmp1{i}(1:rold,:)=[];
0251          core_add=<a href="#_sub1" class="code" title="subfunction [super_core]=compute_supercore(i,elem_fun,d,n,ry,ind_left,ind_right,vec)">compute_supercore</a>(i+1,elem_fun,d,n,rtmp,tmp1,ind_right,vec);
0252          new_core=zeros(ry(i+1),n(i),n(i+1),ry(i+3));
0253          <span class="keyword">if</span> ( rold &gt; 0 ) 
0254            new_core(1:rold,:,:,:)=super_core{i+1};
0255          <span class="keyword">end</span>
0256          new_core(rold+1:<span class="keyword">end</span>,:,:,:)=core_add;
0257          super_core{i+1}=new_core;
0258          <span class="comment">%Recompute supercore{i+1};</span>
0259        <span class="keyword">end</span>
0260        
0261        <span class="comment">%Now one can check, if which supercore is better approximated</span>
0262    
0263        
0264    <span class="keyword">end</span>
0265        <span class="keyword">if</span> ( i ~= 1 &amp;&amp; i ~= d - 1 &amp;&amp; change_dir_on) 
0266          sm=super_core{i-1}; sp=super_core{i+1};
0267           i1m=r_left{i-1}+(i_left{i-1}-1)*ry(i-1); 
0268           i2m=i_right{i}+(r_right{i}-1)*n(i);
0269            i1p=r_left{i+1}+(i_left{i+1}-1)*ry(i+1); 
0270            i2p=i_right{i+2}+(r_right{i+2}-1)*n(i+2);
0271            
0272            sm=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(sm,[ry(i-1)*n(i-1),n(i)*ry(i+1)]);
0273            sp=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(sp,[ry(i+1)*n(i+1),n(i+2)*ry(i+3)]);
0274            smb=sm(:,i2m); [smb,dmp]=<a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(smb,0);
0275            sub_m=smb(i1m,:);
0276            sm_appr=(smb / sub_m) * sm(i1m,:);
0277            
0278            spb=sp(:,i2p); [spb,dmp]=<a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(spb,0);
0279            sub_p=spb(i1p,:);
0280            sp_appr=(spb / sub_p) * sp(i1p,:);
0281 
0282 
0283        <span class="keyword">end</span>
0284    <span class="comment">%Convert iadd1 to i_left and i_right format; compute new</span>
0285    <span class="comment">%ind_left{i+1},ind_right{i+1}</span>
0286    
0287    <span class="keyword">if</span> ( verb &gt;= 1 ) 
0288          fprintf(<span class="string">'Step %d sweep %d ry(%d): %d -&gt; %d rank=%3.1f \n'</span>,i,swp,i+1,ry(i+1)-numel(iadd1),ry(i+1),mean(ry));
0289    <span class="keyword">end</span>
0290    <span class="keyword">if</span> ( i ~= 1 &amp;&amp; i ~= d-1 &amp;&amp; change_dir_on )
0291       er1=<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(sm-sm_appr,<span class="string">'fro'</span>)/<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(sm,<span class="string">'fro'</span>);
0292       er2=<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(sp-sp_appr,<span class="string">'fro'</span>)/<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(sp,<span class="string">'fro'</span>);
0293    <span class="keyword">else</span>
0294      er1=0; er2=0;
0295    <span class="keyword">end</span>
0296        
0297 
0298    <span class="keyword">if</span> ( (strcmp(dir,<span class="string">'lr'</span>) &amp;&amp; er1 &gt; er2) &amp;&amp; er1 &gt; eps &amp;&amp; change_dir_on )
0299         dir=<span class="string">'rl'</span>;
0300         <span class="keyword">if</span> ( verb &gt;= 1 ) 
0301                 fprintf(<span class="string">'i=%d CHANGED DIRECTION LR-&gt;RL er1=%3.2e, er2=%3.2e\n'</span>,i,er1,er2);
0302         <span class="keyword">end</span>
0303 
0304         i=i-1;
0305         rank_increase=true;
0306    <span class="keyword">elseif</span> ( strcmp(dir,<span class="string">'lr'</span>)  )
0307         i=i+1;
0308    <span class="keyword">elseif</span> ( (strcmp(dir,<span class="string">'rl'</span>) &amp;&amp; er2 &gt; er1 ) &amp;&amp; er2 &gt; eps &amp;&amp; change_dir_on)
0309         dir=<span class="string">'lr'</span>;
0310         i=i+1;
0311         <span class="keyword">if</span> ( verb &gt;= 1 ) 
0312           fprintf(<span class="string">'i=%d CHANGED DIRECTION RL-&gt;LR er1=%3.2e, er2=%3.2e\n'</span>,i,er1,er2);
0313         <span class="keyword">end</span>
0314         rank_increase=true;
0315    <span class="keyword">elseif</span> ( strcmp(dir,<span class="string">'rl'</span>) )
0316         i=i-1;
0317    <span class="keyword">end</span>
0318                
0319            <span class="comment">%fprintf('er1=%3.2e, er2=%3.2e \n',er1,er2);</span>
0320    
0321    <span class="keyword">if</span> ( strcmp(dir,<span class="string">'lr'</span>) )
0322      <span class="keyword">if</span> ( i == d-1 )
0323         dir=<span class="string">'rl'</span>;
0324      <span class="keyword">end</span>
0325    <span class="keyword">else</span>
0326      <span class="keyword">if</span> ( i == 1 )
0327         dir =<span class="string">'lr'</span>;
0328         swp=swp+1;
0329         <span class="keyword">if</span> ( rank_increase ) 
0330           rank_increase=false;
0331         <span class="keyword">else</span>
0332           converged=true;
0333         <span class="keyword">end</span>
0334      <span class="keyword">end</span>
0335    <span class="keyword">end</span>
0336 <span class="keyword">end</span>
0337 <span class="keyword">if</span> ( verb &gt;= 1 ) 
0338 fprintf(<span class="string">'Converged in %d sweeps  \n'</span>,swp);
0339 <span class="keyword">end</span>
0340 <span class="comment">%Now compute the approximation itself. Will it be easy?</span>
0341 <span class="comment">%Simple idea: compute cores out of the supercores; they-are-the-same</span>
0342 
0343 
0344 <span class="comment">%% Computation of the approximation</span>
0345 <span class="comment">%Supercore-&gt;submatrix-&gt;SVD (?) -&gt; new core</span>
0346 ps=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0347 cr=zeros(ps(d+1)-1,1);
0348 
0349 <span class="keyword">for</span> i=1:d-1
0350   score=super_core{i};
0351   score=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(score,[ry(i)*n(i),n(i+1)*ry(i+2)]);
0352   i2=i_right{i+1}+(r_right{i+1}-1)*n(i+1);
0353   i1=r_left{i}+(i_left{i}-1)*ry(i); 
0354 
0355   <span class="comment">%Now do just ordinary check</span>
0356   r=0; 
0357   er=2*eps;
0358   <span class="comment">%sbm=score(i1,i2);</span>
0359   <span class="comment">%score=score(:,i2);</span>
0360     <span class="comment">%We have to solve (approximattttely) the equation</span>
0361    <span class="comment">%Newcore*sbm=score</span>
0362    <span class="comment">%sbm can be rank-deficient</span>
0363    
0364   <span class="comment">%Try and test the rank</span>
0365   sbm=score(i1,i2);
0366   [u,s,v]=svd(sbm,<span class="string">'econ'</span>); s=<a href="../../tt2/@qtt_tucker/diag.html" class="code" title="function [qt]=diag(qt)">diag</a>(s);
0367   <span class="comment">%</span>
0368   maxr=max(<a href="../../tt2/@tt_matrix/size.html" class="code" title="function [sz] = size(tt)">size</a>(sbm));
0369   <span class="keyword">while</span> ( r &lt; maxr &amp;&amp; er &gt; eps )
0370      r=r+1;
0371       u0=u(:,1:r); s0=s(1:r); v0=v(:,1:r);
0372      iu=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(u0); 
0373      iv=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(v0);
0374      iu1=i1(iu); iv1=i2(iv); 
0375      w1=score(:,iv1); [w1,dmp]=<a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(w1,0);
0376      w1=w1 / w1(iu1,:);
0377      appr=w1*score(iu1,:);
0378      er=<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(score-appr,<span class="string">'fro'</span>)/<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(score,<span class="string">'fro'</span>);
0379   <span class="keyword">end</span>
0380   <span class="comment">%Now the most delicate part: correctly form the core</span>
0381   
0382   <span class="comment">%cr(pos:pos+ry(i)*n(i)*ry(i+1)-1)=w1(:);</span>
0383   w2=zeros(ry(i)*n(i),ry(i+1));
0384   w2(:,iu)=w1;
0385   cr(ps(i):ps(i+1)-1)=w2(:);
0386   <span class="comment">%norm(w2*score(i1,:)-score,'fro')</span>
0387   <span class="comment">%keyboard;</span>
0388   <span class="comment">%pos=pos+ry(i)*n(i)*ry(i+1);</span>
0389   <span class="comment">%And fix the (i+1)-th supercore</span>
0390   <span class="comment">%iv provides the truncation ry(i+1)-&gt;r</span>
0391   
0392   
0393  
0394 <span class="keyword">end</span>
0395 <span class="comment">%The last one</span>
0396 score=super_core{d-1};
0397  score=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(score,[ry(i)*n(i),n(i+1)*ry(i+2)]);
0398  i1=r_left{d-1}+(i_left{d-1}-1)*ry(d-1); 
0399 score=score(i1,:);
0400 score=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(score,[ry(d),n(d),ry(d+1)]);
0401 cr(ps(d):ps(d+1)-1)=score(:);
0402 
0403 y=<a href="../../tt2/@tt_tensor/tt_tensor.html" class="code" title="function t = tt_tensor(varargin)">tt_tensor</a>;
0404 y.core=cr;
0405 y.ps=ps;
0406 y.r=ry;
0407 y.n=n;
0408 y.d=d;
0409 <span class="comment">%% RANDOMIZED CHECK</span>
0410 <span class="comment">%Now perform a randomized check of the approximation</span>
0411 tmp_fun = @(ind) elem_fun(ind) - y(ind); 
0412 <span class="comment">%s=100;</span>
0413 <span class="comment">%ind=zeros(d,s);</span>
0414 <span class="comment">%for i=1:d</span>
0415 <span class="comment">%  ind(i)=randi([1,n(i)],1,1);</span>
0416 <span class="comment">%end</span>
0417 ind=ones(1,d);
0418 <span class="keyword">for</span> s=1:5
0419   [ind]=<a href="#_sub5" class="code" title="subfunction [ind]=find_fiber_maximum(elem_fun,n,ind)">find_fiber_maximum</a>(tmp_fun,n,ind);
0420 <span class="keyword">end</span>
0421   <span class="keyword">if</span> ( tmp_fun(ind) &gt; 10*eps*elem_fun(ind) ) 
0422     <span class="comment">%This is a restart</span>
0423     restart=true;
0424     <span class="comment">%Enlarge the index set by ind</span>
0425     <span class="comment">%We have to parse the index set</span>
0426     rprev=1;
0427     i=1;
0428     
0429     <span class="keyword">for</span> i=1:d-1
0430        i_left{i}=[i_left{i};ind(i)];
0431        <span class="keyword">if</span> ( i &gt; 1 ) 
0432          r_left{i}=[r_left{i};ry(i)+1];
0433        <span class="keyword">else</span>
0434          r_left{i}=[r_left{i};1];    
0435        <span class="keyword">end</span>
0436     <span class="keyword">end</span>
0437     <span class="keyword">for</span> i=2:d
0438       i_right{i}=[i_right{i};ind(i)];
0439       <span class="keyword">if</span> ( i &lt; d ) 
0440         r_right{i}=[r_right{i}; ry(i+1)+1];
0441       <span class="keyword">else</span>
0442         r_right{i}=[r_right{i}; 1];
0443       <span class="keyword">end</span>
0444     <span class="keyword">end</span>
0445     ry(2:d)=ry(2:d)+1;
0446     <span class="comment">%Bolshie zaglushki</span>
0447     ind_left=<a href="#_sub3" class="code" title="subfunction [ind_left]=get_multi_left(i_left,r_left,ry)">get_multi_left</a>(i_left,r_left,ry);
0448        
0449     ind_right=<a href="#_sub4" class="code" title="subfunction [ind_right]=get_multi_right(i_right,r_right,ry)">get_multi_right</a>(i_right,r_right,ry);
0450     <span class="keyword">for</span> i=1:d-1
0451        super_core{i}=<a href="#_sub1" class="code" title="subfunction [super_core]=compute_supercore(i,elem_fun,d,n,ry,ind_left,ind_right,vec)">compute_supercore</a>(i,elem_fun,d,n,ry,ind_left,ind_right,vec);
0452     <span class="keyword">end</span>
0453     <span class="comment">%keyboard;</span>
0454   <span class="keyword">else</span>
0455      restart=false;
0456   <span class="keyword">end</span>
0457 <span class="keyword">end</span> <span class="comment">%The restart cycle</span>
0458 <span class="keyword">return</span>
0459 <span class="keyword">end</span>
0460 
0461 <a name="_sub1" href="#_subfunctions" class="code">function [super_core]=compute_supercore(i,elem_fun,d,n,ry,ind_left,ind_right,vec)</a>
0462 <span class="comment">%[super_core]=compute_supercore(i,elem_fun,n,ry,ind_left,ind_right)</span>
0463    index_set=zeros(d,ry(i),n(i),n(i+1),ry(i+2));
0464    <span class="keyword">if</span> ( i == 1 )
0465        ileft=zeros(1,0);
0466    <span class="keyword">else</span>
0467       ileft=ind_left{i-1};    
0468    <span class="keyword">end</span>
0469    <span class="keyword">if</span> ( i == d - 1 )
0470       iright = zeros(1,0);
0471    <span class="keyword">else</span>
0472       iright=ind_right{i+2};
0473    <span class="keyword">end</span>
0474    <span class="keyword">for</span> s1=1:ry(i)
0475       <span class="keyword">for</span> s2=1:ry(i+2)
0476         <span class="keyword">for</span> i1=1:n(i)
0477            <span class="keyword">for</span> i2=1:n(i+1)
0478               index_set(:,s1,i1,i2,s2)=[ileft(s1,:),i1,i2,iright(s2,:)];
0479            <span class="keyword">end</span>
0480         <span class="keyword">end</span>
0481       <span class="keyword">end</span>
0482    <span class="keyword">end</span>
0483    M=ry(i)*n(i)*n(i+1)*ry(i+2);
0484    index_set=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(index_set,[d,M]);
0485    <span class="keyword">if</span> ( vec ) 
0486       super_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(elem_fun(index_set),[ry(i),n(i),n(i+1),ry(i+2)]);
0487    <span class="keyword">else</span>
0488        cur_core=zeros(M,1);
0489       <span class="keyword">for</span> k=1:M
0490          cur_core(k)=elem_fun(index_set(:,k));
0491       <span class="keyword">end</span>
0492       super_core=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cur_core,[ry(i),n(i),n(i+1),ry(i+2)]);
0493    <span class="keyword">end</span>
0494    
0495 
0496 <span class="keyword">return</span>
0497 <span class="keyword">end</span>
0498 
0499 <a name="_sub2" href="#_subfunctions" class="code">function [iadd1,iadd2]=enlarge_cross(mat,i1,i2,eps)</a>
0500 <span class="comment">%[iadd1,iadd2,flag]=enlarge_cross(mat,i1,i2,eps)</span>
0501 <span class="comment">%Flag is raised if the initial cross is &quot;very&quot; bad</span>
0502 <span class="comment">%Tests whether the current index set gives a reasonable approximation</span>
0503 <span class="comment">%to the matrix, and enlarges the basis if necessary.</span>
0504 <span class="comment">%The method acts as follows.</span>
0505 
0506 magic_eps=1e-12; <span class="comment">%just some magic</span>
0507 
0508 
0509 <span class="comment">%sbm=mat(i1,i2);</span>
0510 <span class="comment">%[u,s,v]=svd(sbm,'econ');</span>
0511 <span class="comment">%nrm=norm(mat);s=diag(s);</span>
0512 <span class="comment">%r=numel(find(</span>
0513 <span class="comment">%mat_save=mat;</span>
0514 i1save=i1;
0515 i2save=i2;
0516 sbm=mat(i1,i2); 
0517 
0518 <span class="comment">%Compute the initial approximation. The computation is based on the</span>
0519 <span class="comment">%LU decomposition of the sbm until the submatrix is well approximated</span>
0520 [u,s,v]=svd(sbm,<span class="string">'econ'</span>);
0521 s=<a href="../../tt2/@qtt_tucker/diag.html" class="code" title="function [qt]=diag(qt)">diag</a>(s); rm=<a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>(s,min(magic_eps,eps)*<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(s)); 
0522 u=u(:,1:rm); v=v(:,1:rm);
0523 indu=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(u); indv=<a href="../../tt2/core/maxvol2.html" class="code" title="function [ind]=maxvol2(a,ind)">maxvol2</a>(v);
0524 i1=i1(indu);i2=i2(indv); 
0525 
0526 cbas=mat(:,i2);
0527 rbas=mat(i1,:); 
0528 <span class="comment">%rbas=rbas.';</span>
0529 <span class="comment">%[cbas,~]=qr(cbas,0);</span>
0530 <span class="comment">%[rbas,~]=qr(rbas,0);</span>
0531 <span class="comment">%phi=cbas'*mat*rbas;</span>
0532 <span class="comment">%appr=cbas*phi*rbas';</span>
0533 [cbas,~]=<a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(cbas,0);
0534 qsbm=cbas(i1,:);
0535 <span class="comment">%</span>
0536 <span class="comment">% if ( cond(qsbm) &gt; 1e14 )</span>
0537 <span class="comment">%   flag = true;</span>
0538 <span class="comment">%   appr=zeros(size(mat));</span>
0539 <span class="comment">% else</span>
0540 <span class="comment">%   mm=cbas / qsbm;</span>
0541 <span class="comment">%   if ( max(abs(mm(:))&gt; magic_factor ))</span>
0542 <span class="comment">%     flag=true; appr=zeros(size(mat));</span>
0543 <span class="comment">%   else</span>
0544     appr=cbas / qsbm * rbas;    
0545 <span class="comment">%  end</span>
0546 <span class="comment">%end</span>
0547 
0548 
0549 iadd1=[];
0550 iadd2=[];
0551 
0552 <span class="keyword">if</span> ( <a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(mat(:)-appr(:),<span class="string">'fro'</span>) &lt; eps*<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(mat(:),<span class="string">'fro'</span>) )
0553    rnew=numel(i1);
0554 <span class="keyword">else</span>
0555    desirata=eps*<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(mat(:),<span class="string">'fro'</span>);
0556    er=<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(mat(:)-appr(:),<span class="string">'fro'</span>);
0557    mat=mat-appr;
0558    <span class="keyword">while</span> ( er &gt; desirata )
0559       <span class="comment">%Find maximal element in mat</span>
0560       [~,ind]=max(abs(mat(:)));
0561       ind=<a href="../../tt2/core/tt_ind2sub.html" class="code" title="function [ind] = tt_ind2sub(siz,ndx)">tt_ind2sub</a>(<a href="../../tt2/@tt_matrix/size.html" class="code" title="function [sz] = size(tt)">size</a>(mat),ind);
0562       i=ind(1); j=ind(2);
0563       <span class="comment">%if ( ~any(i1==i) )</span>
0564         iadd1=[iadd1;i];
0565       <span class="comment">%end</span>
0566       <span class="comment">%if ( ~any(i2==j) )</span>
0567         iadd2=[iadd2;j];
0568       <span class="comment">%end</span>
0569       u1=mat(:,j); u2=mat(i,:);
0570       u1=u1/mat(i,j);
0571       mat=mat-u1*u2;
0572       er=<a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>(mat(:),<span class="string">'fro'</span>);
0573    <span class="keyword">end</span>
0574 <span class="keyword">end</span>
0575 <span class="comment">% i0=[i1;iadd1]; j0=[i2;iadd2];</span>
0576 <span class="comment">% sbm=mat_save(i0,j0);</span>
0577 <span class="comment">% ss=svd(sbm);</span>
0578 <span class="comment">% fprintf('%10.10e \n',ss)</span>
0579 <span class="comment">% sbm=mat_save(i1,i2);</span>
0580 <span class="comment">% fprintf('AND PREVIOUS \n');</span>
0581 <span class="comment">% ss=svd(sbm);</span>
0582 <span class="comment">% fprintf('%10.10e \n',ss)</span>
0583 
0584 <span class="keyword">return</span>
0585 <span class="keyword">end</span>
0586 
0587 <a name="_sub3" href="#_subfunctions" class="code">function [ind_left]=get_multi_left(i_left,r_left,ry)</a>
0588 <span class="comment">%[ind_left]=get_multi_left(i_left,r_left,ry)</span>
0589 <span class="comment">%Computes (all) left multiindex sets for a given compact representation</span>
0590   d=numel(i_left);
0591   ind_left=cell(d,1);
0592   ind_left{1}=i_left{1};
0593   <span class="keyword">for</span> i=2:d-1
0594       <span class="comment">%ind_cur=zeros(</span>
0595       <span class="comment">%ind_cur is an array of size ry(i-1) x i</span>
0596       ind_cur=zeros(ry(i+1),i);
0597       r_prev=r_left{i};
0598       ind_prev=ind_left{i-1};
0599       i_prev=i_left{i};
0600       <span class="keyword">for</span> s=1:ry(i+1)
0601           <span class="comment">%ind_prev(p1(s</span>
0602          ind_cur(s,:)=[ind_prev(r_prev(s),:),i_prev(s)];
0603       <span class="keyword">end</span>
0604       ind_left{i}=ind_cur;
0605   <span class="keyword">end</span>
0606 <span class="keyword">return</span>
0607 <span class="keyword">end</span>
0608 
0609 <a name="_sub4" href="#_subfunctions" class="code">function [ind_right]=get_multi_right(i_right,r_right,ry)</a>
0610 <span class="comment">%[ind_right]=get_multi_right(i_right,r_right,ry)</span>
0611 <span class="comment">%Computes (all) right multiindex sets for a given compact representation</span>
0612   d=numel(i_right);
0613   ind_right=cell(d,1);
0614   ind_right{d}=i_right{d};
0615   <span class="keyword">for</span> i=d-1:-1:2
0616       <span class="comment">%ind_cur=zeros(</span>
0617       <span class="comment">%ind_cur is an array of size ry(i) x i</span>
0618       ind_cur=zeros(ry(i),d-i+1);
0619       r_prev=r_right{i};
0620       ind_prev=ind_right{i+1};
0621       i_prev=i_right{i};
0622       <span class="keyword">for</span> s=1:ry(i)
0623          <span class="comment">% fprintf('i=%d s=%d \n',i,s);</span>
0624          ind_cur(s,:)=[i_prev(s), ind_prev(r_prev(s),:)];
0625       <span class="keyword">end</span>
0626       ind_right{i}=ind_cur;
0627   <span class="keyword">end</span>
0628   
0629   
0630 <span class="keyword">return</span>
0631 <span class="keyword">end</span>
0632 
0633 <a name="_sub5" href="#_subfunctions" class="code">function [ind]=find_fiber_maximum(elem_fun,n,ind)</a>
0634 <span class="comment">%[ind]=find_fiber_maximum(elem_fun,n,ind)</span>
0635 <span class="comment">%Simple ALS method to compute (approximate) maximum in a tensor</span>
0636 <span class="comment">%In fact, need some non-zero</span>
0637 <span class="comment">%fact=2; %If the new one &lt;= fact times larger than the previous, stop</span>
0638 <span class="comment">%Compute the fibers; find maximum along them;</span>
0639 d=numel(n);
0640 mx=elem_fun(ind); mx=abs(mx);
0641 git=2;
0642 <span class="keyword">for</span> s=1:git
0643 <span class="keyword">for</span> k=1:d
0644   ind_tmp=ind;
0645   <span class="keyword">for</span> i=1:n(k)
0646     ind_tmp(k)=i;
0647     val=abs(elem_fun(ind_tmp));
0648     <span class="keyword">if</span> ( val &gt;= mx ) 
0649        ind(k)=i;
0650        mx=val;
0651        <span class="comment">%fprintf('mx=%f \n',mx);</span>
0652     <span class="keyword">end</span>
0653   <span class="keyword">end</span>
0654 <span class="keyword">end</span>
0655 <span class="keyword">end</span>
0656 <span class="keyword">return</span>
0657 <span class="keyword">end</span></pre></div>
<hr><address>Generated on Wed 08-Feb-2012 18:20:24 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>